Abstract 



In this paper, we present a nonparametric method to estimate the het- 
erogeneity of a random medium from the angular distribution of intensity 
transmitted through a slab of random material. Our approach is based 
on the modeling of forward multiple scattering using Compound Poisson 
Processes on compact Lie groups. The estimation technique is validated 
through numerical simulations based on radiative transfer theory. 
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1 Introduction 

Small-angle scattering of electromagnetic, acoustic or elastic waves constitutes 
the basic tool to characterize porous media pQ, solid suspensions |2j, or Earth's 
lithosphere [3], to cite a few examples only. In many cases, the probing wave is 
assumed to undergo single scattering only and multiple scattering is often con- 
sidered as a nuisance. The main reason is that in the single-scattering regime, 
the diffraction pattern of intensity, i.e. the phase function is related to the power 
spectrum of the medium heterogeneities in a simple manner {4| . In the held of 
optics, some solutions have been proposed to invert for the size of particles in 
the multiple-scattering regime [5] , based on various small-angle approximations 
of the radiative transfer equation [B] . In this work, we develop an alternative ap- 
proach based on the Compound Poisson Process (CPP) model of wave multiple 
scattering in random media. Our model is generic and makes little assumption 
on the underlying wave equation. The random medium is described by its het- 
erogeneity power spectrum which refers to the Fourier transform of the spatial 
two-point correlation function. This enables us to treat on the same footing 
turbid media like Earth's crust or the atmosphere, aggregates of particules, and 
porous media [7] . The finer details of the microstructure encapsulated in higher 
moments of the random field (> 2) are not considered in our approach. 
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The organization of the paper is as follows: In section II, the necessary 
mathematical background on non-commutative harmonic analysis is provided. 
In section III, we develop the Poisson process approach to multiple scattering. 
In section IV, we validate the model through comparisons with Monte Carlo 
simulations of the radiative transfer equation. The inverse problem is examined 
in section V, where we propose an estimator of the power spectrum of hetero- 
geneities and discuss its statistical properties. In section VI, the inverse method 
is validated for realistic examples of random media through numerical simula- 
tions. We consider waves propagating in a turbid medium described by the Von 
Karman correlation function, which is commonly used to represent geophysi- 
cal media. We demonstrate the possibility to estimate the power spectrum of 
the random fluctuations from the angular distribution of intensity transmitted 
through a slab of material. We find that the CPP model is effective before the 
diffusive regime sets in, i.e. when the slab thickness is less than one transport 
mean free path. 

2 Noncommutative Harmonic Analysis 

In this section, we summarize noncommutative harmonic analysis results that 
are important for the present work. We are interested in the harmonic analysis of 
probability density functions (pdf) over compact Lie groups, in order to develop 
a nonparametric approach for the estimation of their characteristic functions 
[8], i.e. their Fourier transform, in Section [5l 

Consider the direction of propagation of a particule or a wave represented 
by the normalized vector \x £ S 2 , where S 2 is the unit sphere in R 3 . We denote 
by /ijv the direction of propagation of the particule after N scattering events. 
Such scattering events are considered random so that fj,N is a random variable 
on S 2 . The relation between the direction after N — 1 scattering events, i.e. 
/Ltjv-i, and the direction after N scattering events /l*jv can be written as: 

/i W = rjvMJV-ij (!) 

where tn is a random element of the rotation group SO (3). This follows from the 
transitive action of 5*0(3) on S 2 . Assuming that a particule enters the random 
medium with initial direction of propagation fiQ, its direction of propagation 
after N scattering events in the medium is given by: 

N 

fi N = r N r N -x . . .r 2 r 1 fi = r„p , (2) 

n=i 

where each r„ represents the random action of the n th scatterer, i.e. the r„ are 
S'0(3)-valued random variables. In the sequel, we present some results about 
the pdf of random variables of the form ([2|) . We make use of known results in 
harmonic analysis on compact Lie groups O I10[ lllj . 
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2.1 Harmonic analysis on SO (3) and S 2 

Let us denote by £ 2 (50(3), dr) the space of square integrable functions on 
50(3) with respect to the Haar measure dr [S]. A probability density function 
(pdf) on 50(3) is a function / e £ 2 (50(3), dr) which obeys the constraint 
f SO t3) f( r )dr = 1. A pdf on 50(3) can be decomposed over the complete 
orthonormal basis of Wigner-D functions Df q (r), with r S 50(3), Z G Z + , 
p,q E Z. This means that given a random variable r S 50(3), its pd/ /(r), can 
be expressed as the infinite series: 

/( r )=E £ (2/+i)/r^rw, (3) 

Z>0 p,q=-l 

where: 

/f 9 = / f(r)Df q (r)dr, (4) 

JSO(3) 

and the overbar means complex conjugation. The set of coefficients ff 9 are 
sometimes called the "Fourier transform" of f(r) [12J . Note that it is possible 
to use a matrix notation, namely f;, in which case, the elements of the (21 + 1) x 
(21 + 1) matrix fj are: (f ; | = ff q . 

I J pq 

Similarly, let us denote by C 2 (S 2 , d[i) the space of square integrable functions 
on the unit sphere S 2 with respect to the invariant measure on the sphere d/i. 
A pdf w € C 2 (S 2 ,d/j,) also satisfies: J g2 w(fi)dfi = 1. Similar to Equation ([3|), 
a pdf w on <S 2 can be decomposed into an infinite series, using the well-known 
spherical harmonics Yf(^x) with /i £ 5 2 , Z € Z + , p € Z. Thus, given a random 
variable /i 6 5 2 , its pdf w(p) can be written as: 

4) = ^^(2i + lgrfM, (5) 

with: 

u)f = / w(n)Y*(n)dp. (6) 

In this case, it is also possible to use a vector representation, where the elements 
of the (21 + 1) vectors W; are wf (p = — Z, . . . , 0, . . . , Z). Again, wf is called the 
Fourier transform of w(fi). 

We now make use of the following fundamental property: the action of 50(3) 
on S 2 is transitive [10] . As a consequence, for any r € 50(3) and /i € 5 2 , we 
have r/i 6 5 2 . This Lie group action of the rotation group on the unit sphere 
implies that, given two pdfs f(r) and w(fi) taking respectively values on 50(3) 
and S 2 , their convolution is a pdf over S 2 and reads: 

h(ji) = (f*w) ( M ) = / f(r)w(r- V)dr, (7) 

JSO(3) 
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for any fi £ S 2 . A very important and interesting consequence of the Peter- Weyl 
theorem [9] is that this convolution equation is transformed into a multiplication 
in the "frequency domain" as follows: 

%)=^^(2i + l)^. (8) 

i>0 p=-l 

Using the matrix/vector notation introduced previously, we obtain the following 
relation: 

hi-f ( W ; VZ, (9) 

which is a matrix- vector multiplication taking place for each degree I. As ex- 
plained in two succesive rotations consist in an other rotation r = r 2 r!. 
Moreover the pdf of r is the convolution of the pdjs of ri and r2, namely /i(ri) 
and f 2(^2), where convolution is again defined with respect to the group action: 

/(r) = / MQMt-^dt. (10) 

JS0(3) 

In equation (|10|). t denotes an element of 5*0(3), and dt stands for the Haar 
measure. Consequently, an N — times product of independent random elements 
of 5*0(3), i.e. N consecutive random rotations, consists in a random rotation 
whose pdf is the N — times convolution of the pdjs of the elementary rotations. 
Thus, if we denote by r = r^r^-i ■ ■ ■ J^2 r i the outcome of N successive rotations, 
the Fourier transform of r can be expressed as follows: 

N 

n vl ( n ) 

n=l 

which is a simple product of matrices. In Ijlip. (^n^j denotes the matrix of 
Fourier coefficients of degree I of f n , the pdf of r„ . 

2.2 Symmetries 

Using Eq. © and (jlip. it is possible to obtain the Fourier coefficients of any 
distribution on S 2 that was multiply convolved by pdfs on SO (3). In what 
follows, we will derive more specific results for distributions on S 2 which are 
functions of fi = cos 8 only, with 8 the scattering angle. Such functions are 
called zonal functions on 5 2 [9]. This simplification relies on the property of 
statistical isotropy of the random medium. As a consequence, the Fourier trans- 
form over 1S 2 reduces to a Legendre expansion over [—1; 1] for zonal functions. 
In anisometric random media where the correlation length depends on space 
direction, this simplification is not permitted. 

In order to obtain simple results, we further impose that the initial distri- 
bution of directions of propagation fj,Q is a zonal function. Such an assumption 
is not too restrictive and covers for example the case of an incident plane wave, 
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i.e. is a Dirac distribution. Using the XZX parametrization of 50(3) [T2"] . 
the symmetry of the scattering process around the direction of propagation 
implies that the pdf w^li) of the direction of propagation /z, parametrized us- 
ing co-elevation 9 (with respect to direction of propagation) and azimuth ip, is 
invariant with respect to <p. This allows us to expand as follows: 

W M = J2( 21 + WiWlti = J2( 21 + l)mPi(cos9), (12) 
;>o z>o 

where the wi are scalar quantities. The rotational symmetry of the scattering 
process also implies that the pdf /(r) of random elements of 50(3) must be 
bi-invariant, i.e. invariant by left and right action on 50(3) (see [3] for details). 
Such a bi-invariant pdf /(r) depends on a single variable \i — cosO, and can 
therefore be expanded as: 

/(r) = £(2/ + l)f?°W¥) = £(2Z + l)M(cos0). (13) 

z>o i>o 

Considering the successive action of N i.i.d. random rotations on the initial 
distribution of propagation directions liq leads to the following expression of 
the pdf Ii(lin) of lin = r N r N -i . . . r 2 ri/i : 

In the special case of an incident plane wave or a highly collimated beam, the 
initial direction obeys the following simple probability distribution w(liq) — 
5(/i — (i z ) where fi z is the vector pointing in the direction of propagation. Intro- 
ducing this expression of w in Equation (|12| yields: wi — Pi(li z ). By symmetry 
ii z must coincide with the co-latitude 0, i.e. the north pole 9 = 0, which implies 
wi = 1, WI. In such a case, eq. (fH|) becomes: 

/(^vH^ + lX/O^M) (15) 

The pdf of /itjv possesses a simple Legendre expansion with coefficients equal to 
the iVth power of the Legendre coefficients fi of the 50(3) random variables. 

Note that Equation [15] could have been obtained using the well-known sum- 
mation formula for Jacobi polynomials (see Chapter 2 in [5]) in the case of zonal 
functions on the sphere. The approach developed in this section is more general 
because it could be applied to less symmetrical distributions, thereby allowing 
the modeling of more complicated scattering processes. 

3 The Compound Poisson Process model for mul- 
tiple scattering 

In this section, we develop the multiple scattering model. We begin with a sum- 
mary of useful results about Compound Poisson Processes (CPP) on compact 
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Lie groups [13] . Note that we use the term pdf when we refer to the function 
describing the angular pattern of scattering. In Section [4} we will substitute it 
with phase function which is the usual terminology in the field of scattering in 
random media. 

3.1 The model 

The implementation of the CPP model for the study of multiple scattering of 
particules/waves has already been described in [T3]. In [2], the authors were 
mostly interested in formulating the direct problem, i.e. in the way a CPP can 
be used to predict the output distribution of scattering angles in a scattering 
experiment. They demonstrated the ability of the CPP to model the multiple 
scattering of electrons and used some recursive integral equations to obtain the 
pdf of the CPP. The CPP model introduced in [T3] is based on the cumulative 
scattering angle, which is a real valued random process defined as a sum of 
real-valued random variables. In this work, we introduce a CPP on a compact 
Lie group, namely the rotation group 50(3). Our approach is more general and 
does not rely on an a priori small-angle approximation. 

Let us consider a particule/ wave entering a random medium at time t = 0, 
with initial direction of propagation /io 6 S 2 , as illustrated in Figure [TJ For 
the moment we neglect the possibility that the wave may escape the random 
medium. We consider a slab of random material with mean free time r and 
velocity c = 1 for simplicity. The normalization of the velocity implies A = 
l/£ where i is the mean free path. We seek to model the time evolution of 
the direction of propagation fi(t). Denoting by N(t) the random number of 
scattering events that have occurred after propagation during time t in the 
random medium, the direction of propagation fi{t) can be written as: 

N(t) 

= II r «A*o- (16) 

n=0 

In order to model wave scattering, N(t) is chosen as a Poisson process, and is 
independent of the r n . To be fully consistent, we must impose ro = I, where I 
denotes the identity in SO (3). This simply means that the unscattered energy 
propagates in the direction imposed by the source. Clearly, the r„, n > 1 are 
independent and identically distributed random rotations described by the pdf 
/(r). Our choice for N(t) implies the usual exponential distribution of times of 
flight between two scattering events, with parameter A = 1/t. 

3.2 Probability density function of the CPP 

To complete our task, we need to relate the pdfoi /j,(t), denoted by p(/j,(t)), to the 
scattering properties of the medium. Keeping in mind the i.i.d. assumption for 
r„ and using conditional probability decomposition with respect to the Poisson 
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process N(t), one gets: 



(\t\ n 

p( f ,(t)) = Y,e~ Xt( -^-f® n (r)*^o), (17) 
' — ' rt! 

n=0 

where ®n denotes the n — times convolution, /(r) is the common distribution of 
the r„ and w(fXo) is the distribution of /xo. Making use of the notation introduced 
in Section |2~T1 for the multiple action of r,: fj, n = r„r„_i . . .r 2 rir /io, the pdf ot 
fi n can be rewritten as: 

/( M „) = /®"(r)* W ( Mo ). (18) 

Assuming again that the distribution of fio is a Dirac at the north pole, the 
distribution of fi(t) simplifies to: 

p(M«))=E^P^*/ 8n (r). (19) 

The assumption about /iq is not very strong since in many practical cases the 
direction of propagation of the incoming wave is either known or controlled. 
For an isotropic random medium, we recall that the pdfs of r; can be consid- 
ered zonal. Therefore, the pdf /(r) can be expanded in a Legendre series with 
coefficients denoted by fk- Using eq. (fT5j) the pdf of /i(t) takes the form: 



pW)) = E^r e " At E^ + 1 )(/0" Pi( ^ ) 



n< 

n=0 l>0 



(20) 



e- At y(22 + l)e A ^P,(^ 



This provides an expression of the distribution of /i(t) of the form: 



P (m) = Y,( 2l + 1 > xtu '~ 1)p ^)- (2i) 



Expanding p(fi(t)) in a Legendre series: 



p(/i(t)) = E( 2Z + 1 )w^(M), (22) 
z>o 

a term-by-term identification in Equations (|2Tj) and ([22]) yields: 

fa = e Xt ^- l \ (23) 

The result (|2"3")) relates directly the Legendre coefficients fii of the pdf of the 
scattering process to the Legendre coefficients /; of the pdf of the random 
rotations caused by the scatterers. This close link between the Fourier coeffi- 
cients of these pdfs is the cornerstone of our approach as it enables us to develop 
a statistical estimator of the Fourier coefficients // , from the observations of the 
pdf of /j,(t). 
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Figure 1: Position of the problem. A plane wave or a narrow beam are normally 
incident on a slab of random heterogeneous material with thickness z expressed 
in transport mean free path units. The transmitted intensity probability distri- 
bution is 1(0). 

4 Validity of the CPP model 

4.1 Models of random media 

In order to show the applicability of our statistical approach, we consider the 
following experiment (see Figure [1]) : A plane wave or narrow beam is normally 
incident on a slab of random material described by Von Karman or Gaussian 
spectra, which are commonly used to describe geophysical media [15] , The 
angular pattern of intensity transmitted through the slab is measured on output. 
In the slab geometry with strong forward scattering, we expect that most of the 
energy will be collected in the forward direction and therefore the CPP should 
give some useful predictions of the pattern of transmitted intensity. In fact, this 
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idea is confirmed by the agreement between the CPP formula (pH) and small- 
angle scattering approximations of the radiative transfer equation [6] . It is to 
be noted that formula pTjl was obtained under the hypothesis that we are able 
to collect all the energy at a given time t but is totally free from a small-angle 
scattering assumption. 

In the CPP approach, the building block for multiple scattering is the pdf 
of the random rotations that scatter the initial direction of the incoming beam. 
In the language of scattering theory, the function that describes the angular 
scattering pattern is called the phase function and will be denoted by /(/i), 
where /i = cos(9) is the cosine of the scattering angle. For sufficiently weak 
fluctuations, the phase function can be obtained from the Born approximation 
[4] and depends of the power spectrum of heterogeneities $ as follows: 

/(/i) = c$(2A;sin(0/2)), (24) 

where k denotes the central wavenumber of the probing wave, and c is a nor- 
malization constant. For the sake of clarity, we review basic results on random 
media of special interest. 

4.1.1 Gaussian Media 

The phase function of a Gaussian random medium is defined by the following 
formula: 

HI \ k 2 d 2 ,22 (l-M) 

= 2 {l-e-^f ' (25) 

where a is the correlation length of heterogeneities. The anisotropy parameter 
or mean cosine of the scattering angle is defined as: 



g = J fMndn. (26) 
-l 

This parameter plays a crucial role in the definition of the transport mean free 
path I* = 1/(1 — g). The transport mean free path is the typical length scale 
beyond which multiple scattering can be described by a diffusion equation. In 
the common case of strong forward scattering it can be much larger than the 
mean free path. As will be shown below, our approach is useful in this regime. 

For Gaussian random media, one obtains g = coth ^-^y - ) — W^i- The Legendre 
coefficients of the phase function are given by: 

fi = J f(riPi(l*W- (27) 

In the high-frequency limit (large ka) , the Legendre coefficients of the Gaussian 
phase function can be approximated as: 

W (i+1)/2 , (28) 

with g w 1 - I /2k 2 a 2 . 
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4.1.2 Von Karman Media 



The Von Karman spectrum implies the following angular dependence of the 
phase function: 

m = 2k 2 a 2 (a-l) 

(l - (4fc 2 a 2 + I) 1 "") (1 - 2Pa 2 ( A i - 1))" ' 

where a is an exponent which controls the small-scale roughness of the medium. 

l 

Note that the phase function is normalized ( J f(fi)dfi = 1). The special cases 

-l 

a = 3/2,2 correspond to the well-known Henyey-Greenstein and exponential 
phase functions, respectively. 

Using integration by parts, the coefficients of the expansion can be expressed 

as: 

l 

fl= 2^! /(!" M 2 )'/ (0 (M)dA*, (30) 

-l 

where /^(/-O denotes the I th derivative of the phase function. Using tables of 
integrals, (|30|) yields the following compact form for the Legendre coefficients of 
the Von Karman correlation function: 

t _ 2^{a - 1)(1 + 2ka 2 )- {l+a) r{l + a) 

rQ + /)r( tt )(i-(i + 4fca a ) 1 - ') 

,'l + a 1 + 1 + a 3 , 4ka 4 

x 2-F1 ( — ^— , — ^ — + 



2 \2 



2 ' 2 '2 (l + 2fca 

where the definition of the hypergeometric function 2-P1 can be found in [16] . 
As noted before, in the case a = 3/2, we recover the Henyey-Greenstein (HG) 
phase function which is a classically used approximation to the Mie theory for 
spherical scatterers. The Legendre coefficients can be put in the form /; = g, 
with g = (1 + 2k 2 a 2 — \fl + 4k 2 a 2 ) /2k 2 a 2 , i.e. the Legendre coefficients are 
simply powers of the anisotropy parameter, valid for any [0, 1[. A remarkable 
property of HG random media is that the convolution on the sphere of unit 
directions of n HG phase functions with same parameter g, is a itself a HG 
phase function with parameter g n : 

00 

/ISO*) =E(2i + l)fl ,n fl(A0- (32) 
1=0 

4.2 Numerical results 

In order to appreciate the strengths and weaknesses of our statistical model 
we show direct comparisons between the CPP and the radiative transfer equa- 
tion for a Henyey-Greenstein random medium in the strong forward-scattering 
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Figure 2: Comparison between Monte Carlo simulations (dots) and CPP model 
for the intensity transmitted through a slab of random medium with Heynyey- 
Greenstein phase function. The four curves correspond to increasing slab thick- 
nes H = 25/8, 25/4, 25/2, 25£, with t the mean free path. The agreement is 
very good in a cone around the forward direction and degrades as the scattering 
angle approaches grazing angles. 
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regime. The anisotropy factor is g — 0.98 and the slab thickness takes the 
values H = 25/8,25/4,25/2,25^?, where £ is the mean free path of the waves 
in the random medium (the transport mean free path £* equals 50£). For sim- 
plicity, we assume that the index of the slab matches that of the surrounding 
medium. The equation of radiative transfer is solved using the Monte Carlo 
code developed by [T7] in the framework of statistically anisotropic media. Fig- 
ure ([2]) shows that the CPP model gives a non-uniform approximation to the 
exact distribution of intensity transmitted through the slab of random material. 
It is noticeable that the CPP model is always wrong in a cone of directions 
perpendicular to the direction of incidence of the wave. The width of the cone 
of direction increases with the slab thickness H . When H becomes of the order 
of the transport mean free path, typically, the CPP model fails to give reliable 
predictions of the angular distribution of intensity. One obvious reason is that 
in the CPP model, we do not prescribe any boundary condition whereas in the 
slab geometry, the intensity must of course vanish in the lower hemisphere of 
propagation directions. 

We further remark that in the slab geometry the distribution of scattering 
events does not follow a simple Poisson law, which also limits the validity of our 
model. To illustrate this point, we examine the validity of the model in Figure[31 
where we show the comparison between the CPP and radiative transfer theory 
for a Gaussian medium with anisotropy parameter g = 0.98. The slab thickness 
takes the values H = 25/16,25/8,25/4,25/2^. The distribution of transmit- 
ted intensity differs markedly from the H-G case (see Figure [2] for comparison). 
Because the Gaussian medium is extremely smooth, the single scattering pat- 
tern is concentrated in a small cone of directions and large-angle scattering is 
excluded. Propagation at large angle can only occur for sufficiently high-order 
multiple-scattering. This entails a significant deviation of the statistics of the 
number of scattering events from the simple Poisson distribution. This point 
is illustrated in Figure [4] where we plot the distribution of scattering events 
around the forward direction (0.8 < \x < 1) and at large angles (0.3 < [i < 0.5) 
for particles propagating through a slab of thickness H — 25/8£ While in the 
forward direction the agreement between the calculated distribution of number 
of scattering events and Poisson statistics is excellent, at large angles a clear 
discrepancy is observed. The distribution is biased towards larger number of 
scattering events and we verified empirically that it is in fact much better fitted 
by a normal distribution. The standard deviation and mean of the Gaussian 
distribution have been adjusted by trial and error. The main point is to illus- 
trate the large deviation from the Poisson distribution. Note however that the 
Poisson law will always be well approximated by a normal distribution (with 
equal mean and variance) if the number of scattering events becomes very large. 

We now discuss briefly the transition towards the diffusive regime. The case 
of strongly anisotropic scattering has been treated analytically in [18] , where it 
is shown that the pattern of transmitted intensity is almost universal, indepen- 
dent of the type of scatterers. In Figure we show the profile of transmitted 
intensity obtained with the radiative tranfer approach with the slab thickness 
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Figure 3: Comparison between Monte Carlo simulations (dots) and CPP model 
for the intensity transmitted through a slab of random medium with gaussian 
phase function. The four curves correspond to increasing slab thicknes H = 
25/16, 25/8, 25/4, 25/2£, with t the mean free path. The agreement is very 
good in a cone around the forward direction and degrades as the scattering 
angle approaches grazing angles. 
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Number of scattering events 

Figure 4: Distribution of the number of scattering events for the intensity trans- 
mitted through a slab of gaussian random medium with thickness H — 3.125£. 
Squares: (0.8 < fx < 1); dots: (0.3 < \i < 0.5). Solid line: Poisson distribution 
with parameter A = 3.125. Dashed line: Gaussian ditribution with mean fi = 7 
and standard deviation a = 2. 

H = £* /2 for the Gaussian and H-G phase functions with the same anisotropy 
parameter g — 0.98. The curves for the two different phase functions are hardly 
distinguishable. For comparison, we have also plotted the simple outcome of the 
diffusion approximation I(fi) — /x(l/2 + 3/x/2). Our numerical study confirms 
the general conclusions of [T5] and also demonstrates that the diffusive regime 
sets in for a slab thickness of the order of £* but not more. In spite of all its 
deficiencies, the CPP model gives a simple, almost analytical solution of the 
multiple scattering problem and explains accurately the angular distribution of 
most of the forward-scattered energy. 

5 Estimation: the inverse problem 

In this section, we consider the problem of estimating the Legendre coefficients 
of the phase function / of the scatterers from the measurement of the output 
angular distribution p(/i(T)) at a given depth in the random medium. Using the 
CPP model, we can construct an estimator for the Legendre coefficients of the 
phase function /;, using a "decompounding" formula. Note that the mean free 
time t or equivalently (thanks to velocity normalization) the mean free path £ 
has to be known in the decompounding approach. Techniques to estimate the 
mean free path of waves in transmission geometry are described for instance in 
[Hi] . The estimators we derive below are extensions of the work done in [201 HI] 
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Figure 5: Transmitted intensity profile through a slab of thickness H = £*/2 and 
anisotropy parameter g = 0.98. Squares: Gaussian random medium. Solid line: 
H-G random medium. Dashed line: diffusion approximation. The difference 
between the Gaussian and H-G medium is only noticeable close to the forward 
direction (/i = 1). 

for real valued random variables, to the specific case of random variables taking 
values over the rotation group and having zonal pdfs. 

5.1 Estimator for the Legendre coefficients 

After proper normalization, the transmitted intensity profile can be interpreted 
as a pdf p{n(T)) , where T is the time required for ballistic waves to propagate 
through the slab. Let us denote by p(/x(T)) the sample pdf which corresponds 
to the data. The sampling resolution of p(p(T)) depends on the acquisition 
system. Typically, a detector has finite aperture and therefore averages the 
intensity over some finite solid angle d£l. This implies that p(fi(T)) is in fact 
discrete and should be more appropriately termed probability mass function with 
the usual normalization 5Z„=i pC/- 4 ™ 1 ^)) = 1 if M samples are available. 

We are interested in estimating the Legendre coefficients of the pdf p(fi(T)) , 
from the observation of p(/i m (T)). This can be achived in two different ways. 
One can cither define directly the estimator on p(/j,(T)), or use p(/i(T)) to gen- 
erate some realizations of fi(T) and define an estimator with these realizations. 
In this paper, we adopt the second approach. To do so, we use a simple linear 
interpolation to evaluate numerically the cumulative distribution function (cdf ) 
and its inverse from the data p(fi(T)). This allows us to generate as many re- 
alizations of jtt(T) as desired, by evaluating the inverse cumulative distribution 
function for a collection of uniformly distributed random numbers in [0, 1]. 
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Let us assume that we have generated K values of n(T) using this technique 
and let us denote them by fJi K (T) with k = 1, . . . , K. The empirical Legendre 
coefficients of p(fi(T)), denoted by fa are given by: 

^ = |X>(AV(T)). (33) 

re=l 

This estimator is unbiased, i.e. E[/x; — fii] = , V I. This means that the mean 
value of our estimator is the real value we are looking for, i.e. fa. The variance 
of the estimator can be expressed as: 



fa - E ifa 



i(E[i?(„)]-£?). (34) 



where E[Pf( M )]=^ 1 p( At (T))Pf(M)dM- 

In the Gaussian case, the variance takes a simple analytical form, in the 
limit case where (g — 1) <C 1. Using a Taylor series expansion of the Legendre 
polynomials near /i = 1 yields Var(fa) = 1 1 " J ^ 1 > (g — l) 2 . For other random 
media, such as Henyey-Greenstein or exponential, an analytical expression of the 
variance is not easily obtained. This is also applies to Gaussian media, away 
from the limiting case mentioned above. Nevertheless, the behaviour of the 
variance can be obtained by numerical integration. In Figure [6l we present the 
variance for three different phase functions: Exponential, Gaussian and Henyey- 
Greenstein. The estimation has been made using K = 1000 samples. It can 
be seen on Figure [HI that the variance of the estimator increases monotonically 
with I for the three different phase functions. While the Gaussian case follows a 
4 th degree polynomial growth, Henyey-Greenstein and Exponential cases follow 
a more or less linear increase with the degree of the Legendre coefficients I. The 
Gaussian case is favorable for low degree coefficients while Henyey-Greenstein 
exhibits the lowest variance for higher degrees. Note that the variance has 
consequences on the ability of the proposed approach to evaluate correctly the 
high-degree Legendre coefficients. This will be illustrated in Section [5721 

Back to the "decompounding" problem, we wish to estimate the Legendre 
coefficients of / from the samples /i K (T) . This is possible by inverting eq. (|23[) 
and using our empirical estimator (|33p . We obtain an estimate of the Legendre 
coefficients of the phase function of the medium: 

7j = ^lnAi + l. (35) 

Equation (|35f is the "decompounding" formula. It implies that the set of Leg- 
endre coefficients fi can be estimated from the set fa which is directly computed 
from the data thanks to eq. (f3"3"| . Note that it is necessary that fa > 0. This 
is verified for small I (/io > m, Vi > 1) but may not be verified for higher de- 
grees, mainly because of the increase of the variance of the estimator fii with / 
(see Figure [S]). As a consequence, it will be necessary to truncate the number 



17 




5 10 15 20 25 



Degree 

Figure 6: Variance of the Legendre coefficients empirical estimator for differ- 
ent phase functions: Gaussian (squares), Exponential (triangles) and Henyey- 
Greenstein (circles) with anisotropy parameter g = 0.995. 

of Legendre coefficients used in the decompounding approach. The truncation 
degree depends on the acquisition set-up and the kind of medium investigated. 
This will be further illustrated in the following Section. It must be emphazised 
that equation ([35)) is central as it allows direct estimation of the heterogeneity 
spectrum of the medium from the output intensity distribution measurement. 

5.2 Inverting for the random medium power spectrum 

As we will now show, the CPP approach constitutes a powerful tool to invert 
the power spectrum of a random medium in the regime: £ < H < £* , partic- 
ularly in the regime £ <ti £*. In order to demonstrate this, we consider a slab 
of material with Von-Karman correlation function. The cases a = 2 (exponen- 
tial correlation) and a = 5/3 (Kolmogorov spectrum) are considered. Synthetic 
data are generated using Monte Carlo simulations of the scalar radiative transfer 
equation in 3D. It is to be noted that the simulations provide the intensity dis- 
tribution averaged over some finite solid angle, as will be the case in a practical 
experiment. In the numerical simulations, wc imposed d£l = 1/200 in order to 
obtain sufficiently smooth results. Due to the limited averaging, the synthetic 
data can be noisy, particularly where the intensity is small but such phenomenon 
is also expected in practice. To obtain the Legendre coefficients of the intensity 
distribution, we process the data as follows. From the numerical simulations, 
we determine the cumulative distribution of the intensity. Although the orig- 
inal distribution can be noisy the cumulative distribution is much smoother. 
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Figure 7: Cumulative distribution function of the intensity transmitted through 
a slab of random material with thickness H = 6.25£ (dashed line) 12. hi (solid 
line) and anisotropy parameter g — 0.98. Top: exponential medium. Bottom: 
"Kolmogorov" medium. 
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Figure 8: Estimated Legendre coefficients of the heterogeneity power spec- 
trum for two slabs of random material with thickness H = 6.25£ (dots), 12. hi 
(squares) and anisotropy parameter g — 0.98. Top: exponential medium. Bot- 
tom: "Kolmogorov" medium. The dashed line shows the exact Legendre coef- 
ficients of the distributions. 
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Examples of empirical cdf for exponential and Kolmogorov spectra are shown 
in Figure [7] From the cumulative distribution, we can evaluate the coefficients 
of the Legendre expansion of the intensity by drawing a sufficiently large set of 
random numbers ran according to: 

ran — icdf(e), (36) 

where e is a uniformly distributed random number and icdf denotes the inten- 
sity inverse cumulative distribution function obtained by linear interpolation 
of the original cdf. Successive application of formula ([55)) and (|3"5")) yields the 

coefficients fa and the desired Legendre coefficients /; of the power spectrum 
of heterogeneities. To obtain accurate estimates we need to draw typically 10 5 
random samples. The comparison between the estimated coefficients and the 
exact Legendre expansion for two different Von Karman random media with 
g = 0.98 is shown in Figure [8] and shows excellent agreement. In order to check 
the accuracy of the decompounding formula, it is important to analyze at least 
two sets of data corresponding to two different slab thicknesses. First, this offers 
a consistency check for the only parameter that enters in the decompounding, 
i.e. the mean free path of the waves. If there is a significant error on this 
quantity (typically more than 20 %), the estimated Legendre coefficients of the 
power spectrum will differ significantly at low degree I. Second, it provides a 
test of validity of the inferred heterogeneity power spectrum at larger degree 

I. Figure [5] shows that the estimated /; for two different slab thicknesses split 
up beyond some degree Iq. For I > l , it is not possible to estimate reliably 
the coefficients of the original distribution. It is to be noted that this test is 
independent of any assumption on the random medium. 

6 Conclusion 

We have developed a non parametric method to infer the properties of random 
media. Our approach relies on a generic mathematical model and could in prin- 
ciple be used to probe random media with acoustic, elastic, or electromagnetic 
waves. The method also provides tests of consistency and accuracy of the re- 
sults. We show that the angular distribution of intensity in a random medium 
can be "decompounded" to estimate the power spectrum of heterogeneities in a 
random medium. An extension of the theory to incorporate polarization mea- 
surements is currently being investigated. 
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